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ABSTRACT 

In  this  paper  we  discuss  the  applications  of  high  order  compact  finite  difference  methods 
for  shock  calculations.  The  main  idea  is  the  definition  of  a  local  mean  which  serves  as  a 
reference  for  introducing  a  local  nonlinear  limiting  to  control  spurious  numerical  oscillations 
while  keeping  the  formal  accuracy  of  the  scheme.  For  scalar  conservation  laws,  the  resulting 
schemes  can  be  proven  total  variation  stable  in  one  space  dimension  and  maximum  norm 
stable  in  multiple  space  dimensions.  Numerical  examples  are  shown  to  verify  accuracy  and 
stability  of  such  schemes  for  problems  containing  shocks.  The  idea  in  this  paper  can  also  be 
applied  to  other  implicit  schemes  such  as  the  continuous  Galerkin  finite  element  methods. 
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1  Introduction 


Compact  schemes  are  methods  where  the  derivatives  are  approximated  not  by  polynomial 
operators  but  by  rational  function  operators  on  the  discrete  solutions.  In  this  paper  we  are 
interested  in  solving  a  hyperbolic  conservation  law 

ut  +  f(u)x  +  g(u)v  =  0 

u(x,2/,0)  =  u°(x,y)  (1.1) 

using  compact  schemes.  In  the  semi-discrete  form,  a  compact  scheme  for  solving  (1.1)  can 
be  written  as 

=  Vfl,s(»))«  =  £(>*)«  (1-2) 

where  A  and  B  are  both  local ,  one  dimensional  operators.  The  subscript  x  or  y  indicates 
that  the  operator  is  applied  in  the  x  or  y  direction. 

For  example,  a  fourth  order  central  compact  scheme  is  given  by  (1.2)  with 

(At>),  =  i(t>i_!  +  4v<  +  v<+i) 

(■ Bv)i  =  i(vi+1  -  t>i-i),  (1.3) 

a  sixth  order  central  compact  scheme  is  given  by 

(At?),-  =  i(u,_r  +  3v,  -|  u,+i) 
o 

( Bv)i  =  ~r(U)+2  +  28u,+i  -  28u,_i  -  u,_2),  (1.4) 

ou 

and  two  third  order  upwind  compact  schemes  are  given  by 

(Av)t-  =  ^(-v,_i  +  5v,  -  u,+i) 

( Bv)i  =  ^(3v,  -  4u<_i  +  u,-2)  (1.5) 

and 

(Av),  =  ~  ut+i) 

( Bv)i  =  ^(-u<+2  +  4ui+1  -  3u,)  (1.6) 

depending  upon  the  wind  direction.  Notice  that  (1.5)  and  (1.6)  have  the  same  implicit  part 
A  which  is  symmetric.  This  fact  will  be  used  later  in  Section  2  to  define  our  local  means. 


The  cost  of  compact  schemes,  regardless  of  the  number  of  space  dimensions,  involves  only 
inversion  of  the  narrowly  banded  (usually  tridiagonal)  matrix  A  and  hence  is  comparable 
to  explicit  methods.  This  is  notably  different  from  other  implicit  methods  such  as  the 
continuous  Galerkin  finite  element  methods  in  multiple  space  dimensions,  even  if  they  are 
similar  in  one  space  dimension. 

The  advantages  of  compact  schemes  include  the  relatively  high  order  of  accuracy  using  a 
compact  stencil  (for  example,  the  fourth  order  scheme  (1.3)  when  discretized  in  time  using 
Euler  forward,  uses  only  a  three  point  stencil  in  each  time  level),  a  better  (linear)  stability, 
and  usually  fewer  boundary  points  to  handle.  In  recent  years  compact  schemes  have  attracted 
much  attention  in  various  fields  such  as  the  direct  numerical  simulations  of  turbulence.  We 
refer  the  readers  to  [2],  [3],  [4],  [12],  and  [18]  for  more  details. 

The  objective  of  this  paper  is  to  apply  compact  schemes  for  shock  calculations.  As 
with  any  other  linear  schemes  (schemes  which  are  linear  when  applied  to  linear  equations), 
compact  schemes  usually  demonstrate  nonlinear  instability  when  applied  to  discontinuous 
data.  We  follow  the  TVD  (total  variation  diminishing)  ideas  in  [9],  [13]  and  try  to  define 
a  suitable  nonlinear  local  limiting  to  avoid  spurious  oscillations  while  keeping  the  formal 
accuracy  of  the  scheme.  Notice  that  the  compact  scheme,  like  an  implicit  scheme,  is  global. 
That  is,  the  approximation  to  f(u)x  at  x  =  aq  involves  Uk  along  the  whole  line  due  to  the 
tridiagonal  inversion  A-1.  Our  main  idea  is  to  define  a  local  mean,  and  use  it  as  a  reference  for 
introducing  a  local  limiting.  In  Section  2  we  introduce  the  limiting  for  one  space  dimension 
and  prove  total  variation  stability.  In  Section  3  we  introduce  the  limiting  for  multiple  space 
dimensions  and  prove  maximum  norm  stability.  In  Section  4  we  present  numerical  examples, 
and  concluding  remarks  are  included  in  Section  5. 

The  ideas  in  this  paper  were  first  used  by  us  for  continuous  Galerkin  finite  element  method 
in  [7].  That  is  an  on-going  project.  In  this  paper  we  restrict  our  attention  to  scalar  problems 
in  order  to  obtain  provable  stability  results.  The  application  of  the  method  to  systems  of 
hyperbolic  conservation  laws  and  to  other  types  of  compact  schemes  (e.g.  [1])  is  currently 
under  investigation. 

In  this  paper,  we  use  the  total  variation  diminishing  (TVD)  Runge-Kutta  type  time  dis¬ 
cretization,  introduced  in  [14],  [17],  to  discretize  the  ODE  in  the  method-of-lines  formulation 
(1.2).  In  the  second  order  case,  the  time  discretization  is 


uW  =  un  +  AtL(un) 

«"+‘  =  +  J“(1)  +  (1.7) 

and  in  the  third  order  case  it  is 

uW  =  un  +  AtL(un) 

u<2>  =  ^un  +  ju(1>  +  ^AtL{uw)  (1.8) 

un+1  =  \un  +  \u(2)  +  ^A<I(u(2)). 

O  J  o 

These  special  Runge-Kutta  type  time  discretizations  are  labelled  TVD  because  it  can  be 
proven  that  under  suitable  restrictions  on  the  time  step  A t  (the  CFL  condition),  the  full 
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discretization  (1.7)  or  (1.8)  is  TVD,  or  stable  under  another  norm  (for  example,  the  L0 o 
norm)  if  the  first  order  Euler  forward  time  discretization  for  (1.2) 


Un+1  =un  +  A  tL{un) 

(1.9) 

is  TVD  or  stable  under  the  other  norm.  For  details,  see  [14]  and  [17]. 

We  thus  only  need  to  consider  the  Euler  forward  scheme  (1.9)  for  stability  analysis  in  the 
subsequent  sections. 

2  One  Space  Dimension 

In  one  space  dimension,  equation  (1.1)  becomes 

ut  +  f{u)x  =  0 

u(x,0)  =  u°(x), 

(2.1) 

the  scheme  (1.2)  is 

(2.2) 

and  the  Euler  forward  time  discretization  (1.9)  becomes 

u"+1  =<  +  A  tL{un)i. 

(2.3) 

Scheme  (2.3)  can  be  easily  written  into  a  conservation  form 

“"+1  =  “?  -  -  K-0 

suitable  for  shock  calculations.  However,  the  numerical  flux  /i"+1 

(2.4) 

is  not  a  local  function  of 

u£  due  to  the  tridiagonal  inversion  A~x.  If  we  define 

Ci 

III 

(2.5) 

then  scheme  (2.3)  can  be  left-multiplied  by  A  to  become 

<+1  «  S?  - 

(2.6) 

When  written  into  a  conservation  form, 

,-.n+l  _  -n  (  in  in  \ 

(2.7) 

this  involves  a  numerical  flux  which  is  a  local  function  of  u£. 

For  example, 

A+*  “  j  (/(tut.)  +  /(«<)) 

(2.8) 

for  the  fourth  order  central  scheme  (1.3), 
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(2.9) 


/j+i  =  ^  (/( t*i«)  +  29 /(«,•+.)  +  29/M  +  /Mi)) 
for  the  sixth  order  central  scheme  (1.4),  and 

/i+l  =  5  (3/(«i) (210) 

and 

/j+}  =  J  (-/Mi)  +  3/MO)  (211) 

for  the  two  third  order  upwind  schemes  (1.5)  and  (1.6),  respectively.  Notice  that  scheme 
(2.7)  resembles  a  cell-averaged  (finite  volume)  scheme  [11].  The  u,  in  (2.5),  like  a  cell  average, 
is  a  local  mean  of  u,  defined  by  Au  in  (1.3)  through  (1.6).  Since  the  computation  of  the  flux 
/|+i  in  (2.7)  involves  the  values  of  u,  a  “reconstruction”  from  u  to  u 

Ui  =  (/4-1u)<  (2.12) 

is  needed.  This  reconstruction  is  global. 

It  is  now  rather  straightforward  to  define  the  limiting.  We  first  write 


with  the  requirement  that 


/(«)  =  f+(u)  +  /  iu) 


(2.13) 


df+(u) 


>0, 


df~(u) 


<0. 


(2.14) 


du  du 

The  purpose  of  this  flux  splitting  is  for  easier  upwinding  at  later  stages.  The  simplest  such 
splitting  is  due  to  Lax- Friedrichs 


/*(«)  =  \  (/(«)  ±«“),  o  =  muax  l/'(«)l  (2-15) 

where  the  maximum  is  taken  over  the  range  of  u°(x).  We  then  write  the  flux  fi+i  in  (2.7) 
also  as 


where  /.*  1  are  obtained  by  putting  superscripts  ±  in  (2.8)  through  (2.11). 
Next  we  define 


(2.16) 


dfh  =  iu  -  /+M;  <*/;.  =  /"Mi)  -  /; 


»+; 


(2.17) 


A  .  * 

Here  df*  i  are  the  differences  between  the  numerical  fluxes  f  j  and  the  first  order,  upwind 
fluxes  /+(u,)  and  /~(u,+i).  These  differences  are  subject  to  limiting  for  nonlinear  stability. 
We  define  the  limiting  by 
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Cr’  =  m(^++i,A+/+(u,),A+/+(ii-.)) 

dfff  =  m(df-i,A+f-(iii),A+r(a,+l))  (2.18) 

where  A+u;  =  v,+i  —  v ,•  is  the  usual  forward  difference  operator,  and  the  (now  standard) 
minmod  function  m  is  defined  by 


m(di 


s  /  smini<i<*  |a«|,  i 
=  \  0, 


if  sign(ai)  =  •  ■  ■  =  sign(ak)  =  s 
otherwise. 


Notice  that  the  limiting  defined  in  (2.18)  is  upwind  biased. 
The  limited  numerical  fluxes  are  then  defined  by 


£(f 1  =  /+(*)  +  Ki'  =  r  (*♦.)  - 


-(to)  _  f~i 


/-(to) 


and 


/(to)  _  f+(m)  ,  ?— (m) 

7,+  i  /i+i  •  ■'.'4-1-  * 


If  we  define  the  total  variation  of  the  mean  u  by 


(2.19) 


(2.20) 

(2.21) 


TV(  S)  =  E  I  (2.22) 

i 

we  have  the  following  proposition. 

Proposition  2.1 

Scheme  (2.7)  with  the  flux  (2.21)  is  TVDM  (total  variation  diminishing  in  the  means) 


TV{un+1)  <  TV(un) 


under  the  CFL  condition 


max  (f+\u)  ~  f  ((u)) 

Proof:  We  follow  Harten  [9]  and  write  the  flux  difference  as 

/$  -  St\  =  -Ci+iA+u.  +  A.jA+fi..., 

where 


= 


£> 


A +/-(«,)  -  +  ^/r(r) 

_ ~  3 _  2 

A+u, 

A+r(«(-i) + ^‘r’  -  rf/i'r1 

_ T  2 _ 1  2 

A+u.-i 


(2.23) 


(2.24) 


(2.25) 


(2.26) 
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The  limiting  in  (2.18)  and  the  properties  of  /^(u)  in  (2.14)  clearly  imply 

Ci+i>0,  D{_k>  0  (2.27) 

and 

At  /r  n  \  ^  At  /-2A+/-(ua)  +  2A+/+(u,) 

S  ^Ci+l  +  D*i>  -  Ai  { - - 

The  last  inequality  is  due  to  the  CFL  condition  (2.24).  TVDM  (2.23)  is  now  immediate 
according  to  Harten  [9]. 

□ 


\  <  1  (2.28) 


In  order  to  obtain  total  variation  stability  for  u,  we  need  the  following  simple  lemma. 
Lemma  2.2 

If  there  are  two  numbers  0  <  6  <  1  and  a  >  0,  which  are  independent  of  N,  such  that 
the  N  x  N  matrix  A  =  (atJ)  satisfies: 

1  N 

!<j<w  ]7~i  -  a’  and  £  Kl  -  6  la«l  ’  7  =  i,  *  •  • ,  N  (2.29) 


(strongly  diagonally  dominance  for  the  transpose  of  A ),  then  the  Lx  norm  of  A~x  is  bounded 
independently  of  N, 


<  — 


(2.30) 


Proof:  Let  A  =  diag(an,  •  •  • ,  a/viv),  B  =  A  —  A  and  C  =  —BA  l.  We  have 


N  N 

IlClk,  =  max  Y'|c,,j  =  max  V] 

*  t 1  i 


<  6. 


Hence  it  follows  that 


ll^lk 


< 

< 


||[(/-C)A]-1|k  =I|A-1(/-C)-1|U1 

IIA-'IU.IK/  -  CJ-'IU,  <  ||A-1||t,r-J— 
a 

1  -  S’ 


□ 
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For  most  compact  methods,  the  matrix  A  satisfies  the  condition  (2.29)  for  Lemma  2.2. 
For  example,  in  the  schemes  defined  by  (1.3),  (1.4),  (1.5)  and  (1.6),  A  satisfies  the  condition 
(2.29)  with  6  =  %,  a  =  6;  6  =  §,  a  =  5;  <5  =  §,  a  =  3  and  6  =  |,  a  =  3,  respectively.  For 
such  compact  schemes,  we  can  now  prove  the  total  variation  stability  for  u. 

Proposition  2.3 

If  a  compact  scheme  (2.7)  satisfies  the  conditions  in  Proposition  2.1  and  Lemma  2.2,  then 
it  is  TVB  (total  variation  bounded).  That  is, 

mO  =  E|“"+  i~<\<C  (2.31) 

t 

for  all  n  >  0  and  At  >  0.  Here  C  is  a  constant  independent  of  n  and  At. 

Proof: 

By  (2.12),  we  have 

t  i 

□ 


This  Proposition  guarantees  convergence  of  at  least  a  subsequence  of  the  numerical  so¬ 
lution. 

We  now  discuss  whether  the  limiting  defined  in  (2.18)  maintains  the  formal  accuracy 
of  the  compact  schemes  in  smooth  regions  of  the  solution.  For  this  we  need  the  following 
assumption. 

Assumption  2.4 

u,  =  ( Au)i  —  Ui  -j-  0(Ax2)  (2.32) 

for  all  u  €  C2. 


□ 

This  Assumption  is  satisfied  by  any  compact  scheme  with  a  symmetric  A ,  for  example 
all  those  listed  in  (1.3)  through  (1.6). 

Under  Assumption  2.4,  it  is  easy  to  verify  by  simple  Taylor  expansions  that 


A+^iiik)  =  /*(«,•)* Ax  +  0(Ax2)  k  =  i  -  l,i,  i  +  1 

=  i/±(«0.Ax  +  O(Axl),  (2.33) 

Hence  in  smooth  regions  away  from  critical  points  (critical  points  are  defined  here  as 
points  for  which  f+{u)x  =  0  or  f~{u)x  =  0),  the  second  and  third  arguments  of  the  minmod 
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functions  in  (2.18)  are  asymptotically  of  the  same  sign  as  the  first  argument  and  half  in 
magnitude.  Hence  the  first  argument  will  be  picked  by  the  minmod  function  (2.19)  for 
sufficiently  small  Ax,  thus  yielding 


(2.34) 

This  guarantees  the  original  high  order  accuracy  of  the  scheme  in  such  smooth,  monotone 
regions.  At  critical  points,  the  accuracy  will  degenerate  to  first  order  as  a  generic  restriction 
of  all  TVD  schemes  (see,  for  example,  [13]).  To  overcome  this  difficulty,  we  use  a  modification 
of  the  minmod  function 


,  *  *  * ,  flfc) 


ai,  i f  |ai|  <  M Ax2 

m(ai,  •  •  • ,  afc),  otherwise 


(2.35) 


where  M  is  a  constant  independent  of  Ax.  This  modification  is  discussed  in  detail  in  [5]  and 

m 

With  this  modification  we  can  obtain  schemes  which  are  formally  of  uniform  high  order 
accuracy  and  equal  the  original  unlimited  scheme  in  smooth  regions  including  local  extrema. 
Moreover,  we  can  prove  the  following  proposition. 


Proposition  2.5 

The  conclusions  of  Proposition  2.1  and  2.3  are  still  valid  for  any  n  and  At  such  that 
0  <  nAf  <  T,  with  TVDM  in  (2.23)  replaced  by  TVBM  (total  variation  bounded  in  the 
means) 


TV{un)  <  C  (2.36) 

where  C  is  independent  of  At,  if  the  minmod  function  m  in  (2.18)  is  replaced  by  the  modified 
minmod  function  m  defined  in  (2.35). 

Proof: 

The  proof  is  similar  to  that  contained  in  [15]  and  [5]  and  is  thus  omitted. 

□ 

The  choice  of  the  constant  M  in  (2.35)  is  related  to  the  second  derivative  of  the  solution 
near  smooth  extrema.  For  details,  see  [5]  and  [15].  The  numerical  result  is  usually  not 
sensitive  to  the  variation  of  M  in  a  large  range. 

In  this  paper  we  only  consider  pure  initial  value  problems.  u°  in  (1.1)  is  assumed  to 
be  either  periodic  or  compactly  supported.  For  initial  boundary  value  problems,  u  in  (2.5) 
is  defined  differently  at  the  boundary,  as  is  the  scheme  (2.6).  The  limiting  (2.18)  can  be 
modified  at  the  boundary  so  that  the  scheme  remains  TVDM  (or  TVBM)  and  TVB  for 
initial  boundary  value  problems.  We  refer  the  readers  to  [5]  and  [16]  for  more  details. 

3  Multiple  Space  Dimensions 

For  notational  simplicity  we  only  consider  the  two  dimensional  case  ( 1 . 1 )-( 1 .2).  Three  space 
dimensions  do  not  pose  additional  conceptional  difficulties.  As  before,  we  only  need  to 
consider  the  Euler  forward  time  discretization 
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We  again  define 


(3.1) 


u 


n+1 


=  u"  +  A  tL(un)ij. 


Uij  =  ( AyAxu)ij  (3.2) 

so  that  scheme  (3.1)  can  be  left-multiplied  by  AyAx  to  become 

fi?1  =  fij  -  ^(A,S,/(U"))0  -  ^(A.BJK)),,.  (3.3) 

Here  and  in  what  follows  we  will  use  the  commutativity  of  Ax,  Ay,  Bx  and  By  so  that  a 
product  can  be  written  in  any  order.  Scheme  (3.3)  can  be  written  into  a  conservation  form 


u 


n+1  =  -  £(/r+Lj  -  £u)  - 


gn .  ,) 


Axv ,+  *’J  Ay 

which  involves  numerical  fluxes  fV\  j  and  g™.  ,  as  local  functions  of  For  example, 

*T  5 » J  T  2  4 


(3.4) 


-  1 

fi+%,j  —  2  Ay  (/(u«+ i,i)  +  f{u>j)) 

ffij+L  =  |-A.(y(Uij+i)+y(tty))  (3-5) 

for  the  fourth  order  central  scheme  (1.3),  with  analogous  definitions  for  the  other  schemes. 
Again,  scheme  (3.4)  resembles  a  cell-averaged  (finite  volume)  scheme  [10].  The  ft,j  defined 
by  (3.2)  is  a  local  mean  of  u,  and  a  “reconstruction”  from  u  to  u 


uij  =  {A-x'A-ylu)ij  (3.6) 

is  needed  to  compute  the  fluxes  fi+ij  and  gtj+i  in  (3.4). 

We  remark  that  the  additional  costs  of  implementing  scheme  (3.4),  comparing  with  the 
original  scheme  (3.1),  are  the  two  local  operators  Ax  and  Ay.  The  major  part  of  the  cost  still 
consists  of  the  two  tridiagonal  inversions. 

The  limiting  to  obtain  nonlinear  stability  can  now  be  defined  in  a  dimension  by  dimension 
fashion;  we  can  use  the  one-dimensional  flux  splitting  (2.13),  for  /(u),  to  write  the  flux 
as 


$ '+2  .i 


(3.7) 


where  /t*  ^  .  are  again  obtained  by  putting  superscripts  ±  in,  for  example  (3.5).  The  remain¬ 
ing  definition  of  the  limiting  parallels  that  in  Section  2,  with  a  dummy  index  j  added  for 
the  reference  y  value.  We  still  start  with  the  differences  between  the  high  order  numerical 
fluxes  and  the  first  order  upwind  fluxes 


df ; 


(3.8) 


and  limit  them  by 
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di;l£  =  ^(dfr^Aif-iu^Air^j))  (3.9) 

where  A+v,j  =  u,+1>j  —  u,j  is  the  forward  difference  operator  in  the  x  direction  and  the 
minmod  function  m  is  defined  by  (2.19).  We  then  obtain  the  limited  numerical  fluxes  by 


££!  =  /+(”.v) + <crj:  Ku  =  rftw  -  df-: rj 


-i”)  _  f-/c. 


(3.10) 


<*•“> 

The  flux  in  the  indirection  is  defined  analogously. 

In  light  of  [8]  this  scheme  cannot  be  TVD  in  two  space  dimensions.  However  we  can 
obtain  maximum  norm  stability  through  the  following  proposition. 

Proposition  3.1 

Scheme  (3.4)  with  the  flux  (3.11)  satisfies  a  maximum  principle  in  the  means 


max  u?.+1  <  max  u" 


'J  - 


(3-12) 


under  the  CFL  condition 


max  (/+'(u))  +  max  (-/  '(«))]  2£~  +  [max  (^+'(u))  +  max  {~9  '(«))]  -  b  (3‘13) 

where  the  maximum  is  taken  in  mintJ  u?  <  u  <  max,,;  u™y 

Proof:  Similar  to  the  development  in  Proposition  2.1,  we  can  write  the  flux  differences  as 


f^m)  _  f{™)  = 

*  (m)  *  (m ) 

9lii~9l/h  = 


Cl+ijA+u,j  + 

—Cij+iAy+Uij  + 


(3-14) 


C,+  ^>0,  Dt_y>  0,  C.J+i>0,  D^_k>  0 

due  to  the  flux  splitting  (3.7),  the  limiting  (3.9),  and 

Ax  (^'+2'}  +  +  Ah  (^'3+5  +  -  1 


(3.15) 


(3.16) 


the  CFL  condition  (3.13). 
We  then  have 
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un+1 

*] 


=  <  +  £  -  A-^'AX_,,)  +  £  (C,v+,  Ai<  -  A,-t A'i?,.,) 


At 


Ax 


Ay 


1  Ax  +  D'-hi)  Ay  (C,*’J+5  + 


“5 


which  implies  the  maximum  principle  (3.12)  because  it?*1  is  written  as  a  convex  combination 
of  u-1.,  u"±xj  and  u"J±1  with  positive  coefficients  which  add  up  to  one. 


□ 

In  order  to  obtain  maximum  norm  stability  for  u,  we  need  a  lemma  similar  to  Lemma  2.2. 

Lemma  3.2 

If  there  are  two  numbers  0  <  6  <  1  and  a  >  0,  which  are  independent  of  N,  such  that 
the  N  x  N  matrix  A  =  (a^)  satisfies: 

1  N 

max  <  a,  and  ]T  |a,_,|  <  <5  ja„| ,  i  =  l,---,N  (3.17) 

3  ^  * 

(strongly  diagonally  dominance  for  A),  then  the  Loo  norm  of  A~l  is  bounded  independent 
of  N 


IM-’lk.  <  YZs- 

Proof:  The  proof  is  similar  to  that  for  Lemma  2.2  and  is  thus  omitted. 


(3.18) 


□ 

For  the  compact  methods  we  consider,  the  matrix  A  is  symmetric.  Hence  the  requirements 
(2.29)  and  (3.17)  are  the  same. 

We  can  now  use  Lemma  3.2  to  obtain  the  maximum  norm  stability  for  u. 

Proposition  3.3 

If  a  compact  scheme  (3.4)  satisfies  the  conditions  in  Proposition  3.1  and  Lemma  3.2  for 
both  Ax  and  A„,  then  it  is  stable  in  the  maximum  norm.  That  is, 

max|u?|<C  (3.19) 

*»«?  * 

for  all  n  >  0  and  A t  >  0.  Here  C  is  a  constant  independent  of  n  and  At. 

Proof: 

By  (3.2),  we  have 
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max 


=  max 


(VVO,.  .|  <  1 1-Ay  1 1  Uoo  1 1 -A*  1  I  Uoo 

-  (rh)VKI' 


u” 

o 


□ 

This  Proposition  does  not  guarantee  convergence,  but  it  at  least  guarantees  that  the 
numerical  solution  will  not  blow  up  due  to  instability. 

Under  the  Assumption  2.4  for  both  Ax  and  Ay,  we  can  again  easily  verify  that  the 
limiting  (3.9)  maintains  formally  the  original  high  order  accuracy  of  the  scheme  in  smooth, 
monotone  regions.  The  degeneracy  of  accuracy  at  critical  points  can  once  again  be  overcome 
by  adopting  the  modified  minmod  function  (2.35)  in  the  limiting  (3.9). 

4  Numerical  Examples 

To  test  the  behavior  of  the  schemes  discussed  in  Sections  2  and  3.  we  use  the  one  and  two 
dimensional  Burgers  equation  with  the  smooth  initial  conditions 


u(x,y,  0)  =  0.3  +  0.7sin(x  +  y).  (4.2) 


Both  are  assumed  to  have  2jr-periodic  boundary  conditions.  The  solutions  will  stay  smooth 
initially,  and  then  develop  shocks  which  move  with  time.  The  exact  solution  to  (4.1)  can  be 
obtained  by  following  the  characteristics  and  solving  the  resulting  nonlinear  equation  using 
Newton  iteration.  The  exact  solution  to  (4.2)  is  that  of  (4.1)  with  x  replaced  by  x  +  y  and 
t  replaced  by  2 1.  These  are  standard  test  problems  for  scalar  nonlinear  conservation  laws 
containing  shocks.  For  comparison  with  finite  difference  ENO  schemes  and  finite  element 
discontinuous  Galerkin  methods,  see  [17],  [5]  and  [6]. 

The  schemes  we  test  are  based  on  the  fourth  order  central  scheme  (1.3)  coupled  with  a 
fourth  order  Runge-Kutta  time  discretization  (henceforth  referred  to  as  the  central  scheme), 
as  well  as  the  third  order  upwind  schemes  (1.5)-(1.6)  coupled  with  the  third  order  TVD 
Runge-Kutta  time  discretization  (1.8)  (henceforth  referred  to  as  the  upwind  scheme).  For 
the  flux  splitting  (2.13)  we  use  the  Lax- Friedrichs  splitting  (2.15).  The  time  step  A tn  is 
taken  to  satisfy  a  CFL  condition 
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(4.3) 


in  one  dimension  and 


Atn 

max  Inf  I  — —  <  0.5 
*  Ax  “ 


max 


A  tn  A  tn\ 
Ax  Ay  J 


<  0.5 


(4.4) 


in  two  dimensions.  When  the  modified  minmod  limiter  (2.35)  is  used,  the  constant  M  is 
taken  as  1. 

We  first  test  the  effect  of  limiters  when  the  solution  is  smooth  but  not  monotone.  In 
Figure  1  we  plot  the  L\  error  versus  number  of  grid  points,  in  a  log-log  scale,  at  t  =  0.6  for 
the  one  dimensional  case  and  at  t  =  0.3  for  the  two  dimensional  case.  In  such  scales,  the  error 
should  be  a  straight  line  with  slope  —  k  for  a  k- th  order  method.  We  can  see  that  the  original 
compact  schemes  and  the  schemes  with  modified  minmod  limiter  (2.35)  (henceforth  referred 
to  as  the  TVB  limiter)  give  the  expected  third  and  fourth  order  accuracy  respectively,  while 
the  schemes  with  the  minmod  limiter  (2.19)  (henceforth  referred  to  as  the  TVD  limiter)  give 
only  second  order  accuracy  due  to  the  degeneracy  at  the  critical  points.  We  can  also  see 
that  both  the  central  and  the  upwind  schemes  work  well  for  this  smooth  problem. 

We  then  test  the  effect  of  limiters  when  the  solution  becomes  discontinuous.  In  Figure  2 
we  show  the  results  of  the  original  compact  schemes  at  t  =  2  for  the  one  dimensional  case. 
We  can  see  over-  and  under-shoots  as  well  as  oscillations,  and  in  this  case  the  result  of  the 
central  scheme  is  much  worse  than  that  of  the  upwind  one.  In  Figures  3  and  4  we  show 
the  results  with  the  TVD  and  the  TVB  limiters.  Apparently  the  limiters  have  stabilized 
the  solution,  as  predicted  by  the  theory.  However  the  result  with  the  central  scheme  is  not 
quite  satisfactory.  In  Figures  5  and  6,  we  show  the  pointwise  errors,  in  a  logarithm  scale,  for 
the  numbers  of  grid  points  N  =  10, 20, 40, 80  and  160.  We  can  see  that  the  central  scheme, 
even  with  the  TVB  limiter,  shows  a  reduced  accuracy  for  quite  a  large  region  around  the 
shock.  This  indicates  that,  for  a  scheme  which  is  globally  oscillatory  (like  the  central  compact 
scheme),  limiters  can  render  it  stable  but  may  kill  accuracy  in  smooth  regions  since  there 
are  oscillations  there  to  suppress.  On  the  other  hand,  the  upwind  compact  scheme  work 
well,  with  bigger  errors  for  the  TVD  limiter  near  the  smooth  extremum  which  is  close  to  the 
shock.  The  errors  for  the  two  dimensional  case  are  similar  and  are  not  displayed.  In  the  last 
plot,  Figure  7,  we  show  the  surface  of  the  two  dimensional  solution  at  t  =  1  with  40  x  40 
points  using  the  third  order  upwind  method  with  TVB  limiting. 


5  Concluding  Remarks 

We  have  discussed  a  general  framework  to  apply  local  limiters  on  compact  schemes  via  the 
definition  of  a  local  mean.  The  resulting  schemes  are  proven  TVB  (total  variation  bounded) 
in  one  dimension  and  maximum  norm  stable  for  multiple  space  dimensions.  Numerical 
examples  show  that  the  base  compact  scheme  should  be  upwind-biased  in  order  to  obtain 
high  order  accuracy  after  limiting  for  shocked  problems. 
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Figure  1:  Lx  error  versus  number  of  grid  points  in  log-log  scale  for  smooth  solutions.  Stars: 
compact  schemes  without  limiter;  squares:  with  TVD  limiter;  diamonds:  with  TVB  limiter. 
1(a):  Third  order  upwind  scheme,  ID 


1(b):  Fourth  order  central  scheme,  ID 
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1(c):  Third  order  upwind  scheme,  2D 


1(d):  Fourth  order  central  scheme,  2D 
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Figure  2:  Compact  schemes  without  limiter  for  shocks.  Pluses:  computed  solution;  solid 
line:  exact  solution. 

2(a):  Third  order  upwind  scheme 


Figure  3:  Third  order  upwind  scheme  with  limiters  for  shocks.  Pluses:  computed  solution 
solid  line:  exact  solution. 

3(a):  With  TVD  limiter 


3(b):  With  TVB  limiter 


i 


Figure  7:  Surface  of  third  order  upwind  compact  scheme  with  TVB  limiter  for  shocks,  with 
40  x  40  points. 
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